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ABSTRACT 



An iteration program is described for the determination of equi- 
librium composition and flame temperature of high energy fuels (HEF) and 
other chemical reactions. 

Tables of thermodynamic functions and equilibrium constants needed 
for the calculations were prepared for sixteen (16) elements and compounds. 
Based on th£ ideal gas relationship, the method of computing values for the 
tables is shown. The CRC 102A digital computer was programmed to compute 
the values of thermodynamic functions and equilibrium constants presented 
in the tables; thereby providing a reasonably consistent set of data. 
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1. Introduction. 



The high energy fuels in present use and undergoing experimentation 
produce high temperatures upon combustion. On the assumption that chemi- 
cal equilibrium exists among the products of combustion, it is possible 
to calculate the composition and the flame temperature, thereby pre- 
dicting the theoretical performance of a propulsion system. For a system 
containing N products of reaction the simultaneous solution of at least 
N+l equations is necessary. These may involve (1) mass balance, (2) 
dissociation and (3) energy and entropy balance. For nominal products 
of complete combustion of stoichiometric quantities the solution is com- 
plex, but at high temperatures, as dissociation occurs and the number of 
products of combustion increases, the calculations become increasingly 
difficult. 

Since equilibrium composition and reaction temperature are interde- 
pendent, a .series of successive approximations must be made beginning with 
a trial temperature. A number of methods for solving these equations have 
been proposed, Ref. (1), (8), (9), using trial and error, matrices, or 
other convergent approximations. Hand solutions are laborious, so it was 
proposed to program the solution on the Computer Research Corporation 102-A 
digital computer at the United States Naval Postgraduate School, Monterey, 
California. This project was completed during fall 1958 and spring 1959 
in partial fulfillment of the requirements for the degree of Master of 
Sciences. The problem was based on the constant pressure combustion of a 
boron fuel containing carbon, hydrogen and boron (C,H,B) at atmospheric 
pressure using air as oxidizer and including up to 16 products of combustion. 
The program is also adaptable to use with other fuels, explosives or any 
chemical reaction that involves the elements Carbon, Hydrogen, Boron, Oxygen, 
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Nitrogen, and Argon. 

Tables containing a consistent set of thermodynamic functions are a 
primary necessity when calculating equilibrium composition and flame 
temperature. For the elements and compounds included in this problem no 
such tables could be found in a single source nor in a readily useable 
form. Hence a large portion of the effort in this project was devoted to 
compilation^of data and calculation of the thermodynamic functions given 
in Tables I - XXI. The sensible enthalpy, i.e. enthalpy based on the 
value at absolute zero temperature (H°-H°); the total entropy S°; and the 
change in free energy of formation ^ F£ , containing enthalpy of formation 
at absolute zero, were tabulated for the temperature range 100° to 9900°R, 
Tables II - XVII. Equilibrium constants were calculated for the 11 reac- 
tions shown and these values tabulated over the temperature range 2000°R 
to ?900 C R, Tables XVIII - XXI. All calculations were based on the ideal 
gas relationship. Values of specific heat at constant pressure C° may be 
easily obtained as first differences in enthalpy. 
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2. Procedure and Discussion 



As a representative high energy fuel (HEF) containing carbon C, 
hydrogen H, and boron B: triethyl - diborane was chosen for use in the 
calculations. When burned with air containing twenty-one percent oxygen 
(21% C^), seventy-eight percent nitrogen (78% and one percent argon 
(1% A), a number of compounds could be produced. Including products of 
di ssociatiqji, 16 elements and compounds, Table I, were chosen to work 
with. Based on methods in Refs. (1), (2) and (3), combinations of these 
products of combustion produced 11 possible equilibrium reactions, with 
some repetition shown by the water-gas reaction. 

A readily available source of tabulated thermodynamic functions did 
not exist so it was decided that a reasonably consistent set of tables of 
data would be produced. Formula mass in Table I represents the latest 
values from the International Standards adopted in 1956. The enthalpy of 
formation at absolute zero temperature (H^ ) 0 , Table I, was taken from 
Refs (4) and (6). A program was written for the CRC -102A digital com- 
puter such that basic values of thermodynamic properties could be stored 
in the computer memory and these used to compute the properties at inter- 
mediate temperatures. The enthalpy function was computed for a series of 
small temperature increments by the following formula from Ref. (1): 

H = H fc • C p (T-T b ) 

where the subscript b refers to a value at an arbitrarily selected basic 
temperature near the temperature of concern. Since enthalpy was based on 
the value at absolute zero temperature, the tabulated function is a 
"sensible enthalpy" and is a standard value at any given temperature minus 
the standard value at absolute zero, H*-Hg. 
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The entropy was computed in much the same manner by the formula, 

Ref (1) 

L 

S = S. + C In T. 
b P b 

but since the entropy at absolute zero is equal to zero, the tabulated 
values become the standard value of total entropy S°, at any given tem- 
perature. As will be shown later in the flame temperature calculations, 
the fact that entropy equals zero at absolute zero makes this a good ref- 
erence temperature to use during these calculations. 

Basic values for Tables II - XVII were found in various dimensionless 
and dimensional forms in Refs. (4) - (9) which represented the most con- 
sistent sources available. In the computer these values were all con- 
verted to the English Standard of Engineering units (BTU/#mol) and all 
computations made and results presented in these units. 

Enthalpy H, entropy S, and specific heat at constant pressure for 

A, CO, CO 2 , H. N, N 2 , 0 and O 2 were taken from Refs. (4) and (5): for 

B. BO, 02^3' fi nd from Ref. (6): and C graphite) from Ref. (7). 
Basic values were taken from the references at temperatures of 100°K, 

2°8. 16°K, 400°K, 500°K, 700°K, 1000°K, 1200°K, 1500°K', 2000°K, 2500°K, 
3000°K, 3500°K, 4000°K, 4500°K and 5000°K. Smaller increments were chosen 
in the low temperature range because that is the region in which changes 
more rapidly. A graphical extrapolation was made for Carbon to obtain 
basic values at 4500°K and 5000°K. 

The specific heat Cp may be obtained by taking the first difference 
in enthalpies, i.e. for CO, Table VII: 

H = H 5500°R ~ H 5400°R ~ C p 

= (44,834 - 43,944) BTU/#mol = C p (5500-5400) 

891 BTU/#mol 

C = 100 °R -8.91 BTU/(#mol °R) 

P 
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The change in free energy of formation was calculated on the basis of 
formation from the elements and includes the enthalpy of formation at 
absolute zero temperature, Table I, hence tabulated values are standard 
values (^FJ) for any given temperature- From Chapters 7 and 12 Ref. 

( 1 ) 

F f = H - TS 

* A F f = A H - A (TS) 

but since calculations were made on the basis of enthalpy and entropy 
functions at a single temperature, the term on the right of the equation 
becomes TAS and 

A F° = AH 0 - TAS 0 

The function A F is the Gibbs’ function or Gibbs' free energy, named 
after Willard Gibbs who first established this useful relationship. As 
an illustration; for the formation of atomic nitrogen from the free ele- 
mental molecular nitrogen 

1/2 N 2 = N 
AF° = AH 0 - TAS 0 

AF° = (H£ - 1/2 H£ ) + (H°), N - T(S° - 1/2 S° ) 

‘2 ‘ — 

Tabulated values of enthalpy, entropy, and change in free energy of form- 
ation are presented in Tables II - XVII. 

The equilibrium reactions in Tables XVIII - XXI were formulated and 
equilibrium constants calculated on the basis of methods found in Refs. 
(1). (2) and (3). An important use of the Gibbs* free energy function be- 
came apparent in these calculations. 

A F° = - RT In K 
f P 

which states that when reactants in their standard states are changed to 
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products in their standard states there is a change in free energy equal 
to -RT times the natural logarithm of the equilibrium constant, Ref. (2). 



for convenience in practical problems, the actual values of equilibrium 
constants were tabulated rather than a logarithm function. This limited 
tabulation in the low temperature range, but, with the free energy functions 
and the above formula, the equilibrium constant nay be easily determined 
for any temperature. 

Values in Tables II - XXI were computed on the digital computer and 
punched on IBM cards which were subsequently used in an IBM-402 Accounting 
Machine to print the tabular format. 

With reasonably consistent values of the thermodynamic functions, the 
problem of equilibrium composition and flame temperature was attacked. 

Since these two parameters are interdependent, a trial temperature must be 
chosen and the composition determined; then a check for energy, entropy 
balance must be made and if the temperature was not correct a new trial 
value must be assumed. This procedure must be continued until the correct 
temperature, hence proper composition, is found. 

For the 16 products of combustion, simultaneous solution of 17 equations 
is necessary. An iterative scheme of solution is used on the digital com- 
puter and this is best shown by an illustrative example. First choose a 
trial temperature, 5400°R, and having values of K , Tables XVIII - XXI, the 
interrelation between combustion products and reactants is established and 



Then 





equilibrium components may be determined. Burn the fuel completely with 



stoichiometric quantity of air as oxidizer to determine the amounts of 
nominal (major) products. With triethyl-diborane the following reaction 
occurs 



B 2 H 3 (C 2 H 5 } 3 + 12 °2 + ( 12 > ff N 2 + (12) jTf A = 6 C0 2 + 9 H 2 0 

» + B 2 0 3 + 44.6 N 2 + 0.571 A 

Now since complete burning rarely occurs in an engineering application 
without an excess of oxygen, include CO, H 2 and 0^ among the nominal pro- 
ducts. The reaction then takes the form 



B o H ,(°JO, + 12 0 o + (12) N 0 + (12) A = a C0 2 + b H 2 0 



2 3 2 5 3 



Zl 



2 ! 



, + c B 9 0, +dN+eA+fC0+gH+h0 

z j 2 2 2 

where the letters a-h represent the unknown molar amounts. By a material 
balance equations may be written to equal the number of elements involved. 



a 

b 

c 

d 

e 

f 

g 

h 



CO, 



H 2° 

B 2°3 

N_ 



A 

CO 

H 2 

0„ 



Carbon (C) 
Hydrogen (H) 
Boron (B) 
Nitrogen (N) 
Argon (A) 



a + f = 6 
2b + 2g = 18 
2c = 2 
2d = 44.6 
e = 0.571 



Oxygen (0) 2a + b + 3c + f + 2h = 24 

Since only six (6) equations exist for these eight (8) unknowns one must 
resort to equilibrium equations for solution of the other two. To solve 
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for f choose the following equilibrium reaction from Table XVIII 



C0 2 ~± CO + 1/2 0 2 

then by partial pressure, p, notation 




% 

The partial pressure of a gas may be represented by the mole fraction 
times the total pressure. The unknown number of moles of carbon monoxide 
has been assigned the value f; let, sigma, , equal the total number of 
moles of combustion products, then the mole fraction of CO is f/^. • 

With the value of at 5400°R from Table XVIII the above equation becomes 



K = 0.32898 = 
P 



Vi P( h MP) 



V 



a. 



!z P 



0.32898 




From the carbon (C) material balance, a = 6 - f and 



0.52.1 



i pvt pf 

6 - f 



r __6jo^zT}__ r yj± 

f (VS Pt*i-0.3Z1 (Vs P) h +0.5^ 



where P is the total pressure in atmospheres. Solution for g, the moles 
of H 2> is completed in a similar manner. Then for each new product of 
combustion a new equilibrium equation must be solved. This points out 
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the need for values of equilibrium constants at various temperatures. 

It should be noted that the solution for f contains the unknown 
h, for the number of moles of unburned oxygen; this will also occur in 
the solution for g. Since simultaneous solution of this large number of 
equations is difficult these are best solved by an iterative process. A 
trial value of h is chosen; solution for the other unknowns is made on 
the basis of this trial, then by the summation shown for the oxygen balance 
the trial value is checked. Considerable assistance is gained by making 
a table to show the operations involved, (see Iteration Table page 10). 




Carbon for this example. For products of dissociation involving other 
minor components, assign the following unknowns. 

j = OH, k = H, 1=0, m = N, n * NO, q = BO. 

For the first iteration, only major products are considered to exist 
hence all minor products will be zero. The role played by pressure is 
indicated in the Iteration Table, but for convenience was taken as one (1) 
atmosphere. 

At the end of the first iteration the value of h is determined and 
compared with the trial value. The discrepancy in h is taken as the cal- 
culated value minus the trial value. For the second trial value of oxygen 
divide the discrepancy by two (2) and add this to the first trial value 
of h. Now, however, all minor products of combustion should be included. 

After the second iteration extrapolate (or interpolate) the dis- 
crepancy to zero to aid in the estimate of a trial value of h for the third 
iteration. The third trial should then give a reasonable answer provided 
no errors in arithmetic were made. Once programmed, this sort of iterative 
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ITERATION TABLE 



h 


* trial 0 

? 


2.50 


2.46 


2.50 




C 

* 33.37+h+j+k+l+ra+n+q 


41.37 


41.33 


46.704 




(h/£ P) V * 


.246 


.244 


.231 


f 


1.974 

(h/£ py /2 + 0.329 

% 


3.440 


3.450 


3.530 


a 


- 6 - f 


2.560 


2.550 


2.470 


g 


9 

21.96(h/£ P)**+ 1 


1.410 


1.420 


1.484 


b 


■9-g-j-k 


7.590 


7.580 


4.336 


c 


* 1 + q 


1.000 


1.000 


1.000 


d 


= 22.3 - ra - n 


22.330 


23.300 


21.337 


e 


® 0. 571 (const. ) 


0.571 


0.571 


0.571 




(g/* P) ,/A 




0.185 


0.178 


J 


0.0476 b 

(g/Z P)"*- 


0 


1.960 


1.160 


k 


3 ^/P (0.0248) (g/£ P) ,/2 - 


0 


1.200 


1.310 


1 


3 Z/P (0.0144) (h/£ P)^ 


0 


1.210 


1.300 


m 


- £/p (0.00137) (d/Z pf* 


0 


0.043 


0.043 


n 


fc> J/o 

= 0.1216 d z h ** 


0 


0.920 


0.890 


q 


(c/£ P)^ 

- 2 /p ( 0 . 000039 ) p ;, v 


0 


0.001 


0.001 


h 


2a+b+3c+f+j+l+n+q 
* 12 - « 2 > 


2.42 


0.39 


2.43 


Discrepancy in h s calculated 

- trial 


-0.0s 


-2.07 


-0.07 
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process is quickly and accurately completed on a digital computer. 

Once the composition in number of moles of each product of com- 
bustion has been determined, the trial temperature must be checked. 

This is accomplished by an energy, entropy balance. If the temperature 
is incorrect a new value must be chosen and a new composition determined. 
It is advisable to choose the second temperature value such that it will 
be on the opposite side of the true value from the first trial. Then 
a more accurate estimate may be made and the theoretical flame tempera- 
ture determined. 
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TABLES OF 



THERMODYNAMIC PROPERTIES AND EQUILIBRIUM CONSTANTS 

Symbols 

F Gibbs’ free energy of formation, BTU/#Mol 
H Total heat content, Enthalpy, BTU/#Mol 
Enthalpy of formation, BTU/#Mol 
Kp Equilibrium constant, Dimensionless 
M Formula mass, #/#Mol 

R Gas constant, 1.9872 BTU/(#Mol °R) 

S Entropy, BTU/(#Mol °R) 

T Absolute temperature, °R 

Superscript, denotes standard value at given temperature, 
i.e., H° is standard enthalpy 

o Subscript, denotes temperature, i<>e., H| is standard enthalpy 
at 0°R 



A-l 



TABLE I 



FORMULA MASS AND ENTHALPY OF FORMATION 
AT ABSOLUTE ZERO TEMPERATURE 



Element or 



Compound 


M 


( H f)o 


A 


39.944 


0 


B 


10.820 


0 


* 






BO 


26.820 


- 10,800 


B 2°3 


69.640 


- 543,600 


c 


12.011 


0 


CO 


28.011 


- 48,963 


C °2 


44.011 


- 169,143 


H 


1.008 


92,916 


H 2 


2.016 


0 


h 2 o ( g ) 


18.016 


- 102,793 


N 


14.008 


153,216 


N 2 


28.016 


0 


NO 


30.008 


38,659 


0 


16.000 


105,455 


°2 


32.000 


0 


OH 


17.008 


18,000 



A-2 



TABLE II 



ARGON (A) 



T 


°R 




H°-HS 




s 
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A 

f 


1 


o o 
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-4 
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o 
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